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Abstract 

Runge-Kutta schemes have been used as a method of solving the Euler 
equations exterior to an airfoil. In the past this has been coupled with 
central differences and an artificial viscosity in space. In this study we 
couple the Runge-Kutta time-stepping scheme with an upwinded space approxima- 
tion based on flux**vector splitting. Several acceleration techniques are also 
considered including a local time step, residual smoothing and imiltigrid. 
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Introduction 

We wish to solve the steady-state multidimensional Euler equations with a 
method that Is suitable for a large range of Mach numbers. Al the same time 
we wish the method to be accurate, robust and capture shocks without excessive 
smearing. We also wish to reach the steady state rapidly. To achieve these 
goals we combine the Runge-Kutta scheme Introduced in [1] with the flux-vector 
splitting Introduced in [2]. 

Consider the two-dimensional system 



w^ + f + g * 0. 
t X *y 


We advance the numerical solution In tiiLc using a N stage algorithm 




where f and Dy g are difference approximations to the flux derivatives. 
To check the stability we freeze coefficients and Fourier transform. The 
amplification matrix of (2) la then 





G = 1 +8^ z B 2 + •*• + 8j^ z*^ 

6^ = 1 (3) 

®k “ ^k-1 “N-k+1 

where z is the Fourier transform of At[D f + D gj. When central 

X y 

differences are used then z lies on the imaginary axis. With upwind 

differences z lies on some curve in the negative real half of the complex 
plane. 

Experience has shown that one should usually choose the parameters so that 
the time step is maximal. For central differences this implies that 

CAt/Ax < N - 1 (4) 

where C depends on With upwind schemes no general rules have been 

developed thus far for optimal parameters. At present the parameters have 

been chosen by experimentation. One possibility is presented in the result 

section. 

To appreciate the connection between central differences and flux vector 
splitting we consider a one-dimensional example. Assume that we wish to 
compute a numerical flux at the cell interface i *f V 2 • Quadratlcally 
interpolating the cell averages In zones 1-1, 1, i+1 yields a left-side 

estimate of the state variables 
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(wi^l - 2w^ + 


Interpolating the cell averages in zones i, i+1 , 1+2 yields the right-sided 


estimate 


”i+V2~ ”i+l ~ ^ ^”i+2 " '*■ 12 ^”l+2 ‘ ^”i+l 


The difference between these values Is 0((Ax) J# We now introduce an upwind 
bias in the numerical flux by using from the components of w only those 
characteristic combinations that are converged forward and from 
convected backward. In the approximation of flux splitting this becomes 


rf L ^ \ 




This can be rewritten as 


VV2^ ‘I ^^1+1/2 ^i+V2^ “i‘^i+l/2^VV2“VV2^’ 


dt'*' _ ^ . 
dw dw ' 


L R 

Ql+ V 2 ** Roe-type (3) average of Q over the interval 1/2 » '' 1 + 1 / 2 )* 

From (8) we see that the upwind-biased flux deviates from the average 




flux, used for central differencing, by a third-order term. This leads to a 
fourth-order viscosity with a matrix-valued coefficient. This viscosity 
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prevents the checkerboard Instability similar to the fourth-order viscosity 
introduced In (!]• 

Equations (5) and (6) are modified before their actual use. The first 
order term is multiplied by a switch described in [4] while the second-order 
term is multiplied by its square. When is large compared with 

and ” '*1+1' ^ shock profile, then the limiting yields 

VV2""i+V2“^i+l ”"!• 

Therefore, the viscosity term in (8) now leads locally to a second-order 
viscosity which guarantees a monotone profile. This is similar to the second- 
order artificial viscosity introduced in [1] and discussed in more detail in 
[5]. 

Limiting the high-order terms combined with upwind differencing is a 
robust way of preventing numerical oscillations near discontinuities. To 
achieve the same effect the viscosity of [1] would have to be raised to the 
level of the spectral radius of Q leading to excessive smearing. In 
practice using the code of [1] the opposite approach is used. The coefficient 
of viscosity is adjusted, by trial and error, so that the shock profiles are 
sharp and spurious entropy production in the smooth flow is minimized. For 
smooth flows one can achieve viscosity levels that are smaller than that of 
the upwind scheme with the limiter* For violent flows the artificial 
viscosity is much too large. In any case the coefficient of viscosity used in 
II] is very problem-dependent while the present code has no adjustable 
parameters that need to be played with. 
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Flux splitting In an arbitrary multidimensional body-fitted coordinate 
system can be reduced to a one-dlmenslonal problem with reference direction 
normal to the cell face where fluxes are computed. The resultant code Is a 
full two-dimensional code and does not use time splitting to combine the 
coordinate directions. In all the Interpolations we have Ignored the 
geometric variation of the cells as was done In [1]. This can lead to errors 
for high angles of attack, especially near the trailing edge where cell shape 
and size vary strongly. 

Having discussed Che spatial discretization we now Introduce several 
techniques to accelerate the convergence Co a steady state. The first 
technique Is to use a local time step. This Improves Che running time by an 
order of magnitude due to Che small cells near Che airfoil. The second 
technique Is to use residual smoothing after each stage of the Runge-Kutta 
method. This was first Introduced in [6] for the Lax-Wendroff scheme. If one 
uses central differencing Chen the scheme Is unconditionally stable when the 
smoothing Is done after every even stage. Using an upwind scheme the 
smoothing should be done after each stage of Che algorithm. Even though the 
resultant scheme Is unconditionally stable It Is not efficient to use time 
steps that are Coo large* Time steps about three times as large as those of 
the explicit scheme seem to be optimal. Since Che residual smoothing adds 
only about lOZ to Che cunning time per time step the use of the residual 
smoothing Is advantageous. 

A third acceleration technique Is to use a multigrid method. Jameson [7] 
has proposed using the Runge-Kutta scheme coupled with central differences as 
a smoothing algorithm for a multigrid scheme. The parameters are now 



chosen to damp the high frequencies rather than achieving a maximal time step. 
We use the same technique with upwind differencing in the Runge-Kutta scheme 
instead of central differences with an artificial viscosity. In the central- 
difference version [7] the artificial viscosity is increased compared with the 
standard Runge-Kutta scheme [1]. In addition the raultigrid central-difference 
code seems to rely on the enthalpy damping in order to achieve rapid 
convergence to the steady-state. With the upwind code there is no longer an 
artificial viscosity that can be tuned to damp the high modes. Furthermore, 
the steady-state total enthalpy is not preserved by the flux-vector-splitting 
scheme [2]. Hence, the enthalpy damping introduced in [IJ and [5] cannot be 
used. Nevertheless the multigrid scheme does work with the upwind-biased 
scheme. Thus, the viscosity that is Implicit In the scheme seems to be 
sufficient to compensate for both the artificial viscosity and the enthalpy 
damping of the central-difference scheme. However, the convergence rate of 
the upwind scheme is slower than that of the central-difference scheme, mainly 
because of the enthalpy damping In the latter. 


Results 

The upwind biased version of the code has been run on several different 
cases. The first case Is a NACA-0012 airfoil with M « 0.8, and a » 1.25^. 
The mesh Is an 0-mesh generated by a sheared parabolic transformation. We use 
a four stage Runge-Kutta scheme with Oj - .17, Oj “ *273, - .5 and 

> 1. The residual smoothing is applied after each stage with <• .9X/8 
and B ■ .63 (X is the local Courant number) with X ■ 3. Using a 64 * 16 

y * 




mesh Che residual Is reduced by 4 orders of magnitude after 600 steps. The 
Cp curve Is similar to that achieved by the central-difference code [ 1 ] 
except chat the shock profile is now sharper with one point In the middle of 
the shock along Che airfoil. We have also run several supersonic flows about 
the NACA-0012. The upwind version of the code converges for a larger range of 
Mach numbers Chan does the central difference version. 

Because of Che flux splitting and upwind logic entering Che execution Che 
upwind code Is about two times slower per time step chan Che central- 
difference version. A further slowdown Is caused by the Runge-Kutta method 
which seems to favor the spectral distribution of central differences and 
which has not yet been optimized for upwind differences. Hence, the time step 
Is about half of that for central differences. In addition, Che enthalpy 
damping described in [1] and [3] cannot be used. Thus, Che present version of 
Che upwind scheme Is about 5 times slower in reaching the steady state than 
the central-difference code of [ 1 ]. 

The nulclgrld version of Che code has also been run using a four-stage 
Runge-Kutta. The original parameters were reasonable but a better set is 

» .15, 02 ■ .3275, O 3 “ *57 and o^ » 1. We have also used a six-stage 
formula with o^ ■ .073, 02 * .138, o.^ ■ .22, » .334, Oj ■ .5 and o^ « 1 . 

On a 64 X 16 mesh the nulclgrld version requires fewer iterations to converge. 
However, accounting for the extra work of the mulcigrid Che two codes have 
about Che same convergence rate per work unit. Nevertheless, if a coarser 
mesh is used to initialize Che finer mesh, Chen the total number of supersonic 
points is predicted within 60 Iterations on Che fine mesh. Moreover, it is 
expected that for finer grids Che mltigrid will be more efficient. The 
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central-dlfference nultigrid code is Che fastest code. However, this relies 
heavily on Che enthalpy damping. Thus, for the Navier-*Stokes equations it 
cannot always be used. In this case the upwind multigrid scheme will be the 
most efficient. 
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